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INTRODUCTION 


Recent advances in both computer hardware and numerical techniques have 
led to a significant broadening of the practical choices available for 
analyzing a wide variety of viscous flow problems. Prior to recent 
computational advances, most predictive techniques were confined to inviscid 
flow analyses possibly combined with boundary layer corrections. More 
recently efforts have focused upon more complex flows which were not 
necessarily suited to boundary layer type analyses. One important problem of 
this type is the two-dimensional isolated airfoil flow field problem. 

The two-dimensional isolated airfoil flow field presents a problem which 
has long been of practical interest. The accurate knowledge of airfoil lift, 
drag and moment coefficients under a range of steady and unsteady operating 
conditions is required in assessing the airfoil performance. In the general 
case the airfoil flow field presents complex phenomena even when three 
dimensional effects are neglected. This problem contains viscous regions 
which are laminar, transitional and turbulent, may exhibit extremely strong 
favorable and adverse pressure gradients, may contain multiple regions of 
large separation as well as shed vortices and exhibit important unsteady flow 
characteristics. A particularly complex airfoil flow field occurs in the 
helicopter rotor which may periodically undergo dynamic stall. Dynamic stall 
differs from its static counterpart in two major ways. First of all, the 
maximum lift obtainable under dynamic conditions is greater than that under 
static conditions. Secondly, though under static conditions lift and moment 
are uniquely related to incidence, under dynamic stall conditions the flow 
depends upon the time history of motion and the lift and moment coefficients 
have hysteresis loops associated with them. As the helicopter blade travels 
through the rotor disc, the blade experiences a varying incidence angle. 

Over most of the disc the blade will be unstalled (i.e., the flow will not 
contain any large separated regions leading to a decrease in blade lift or a 
generation of large blade moment coefficients); however, over a portion of 
the disc large regions of separated flow may appear and over these regions 
the blade performance will deteriorate. This time history dependence of the 
problem makes dynamic stall prediction a particularly difficult task. 

Typical experimental investigations of this complex phenomenon are discussed 
in Refs. 1-8. 
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In regard to predictive procedures for airfoils in both steady and 
unsteady flow, airfoils at low incidence are expected to have small viscous 
displacement effects. In these cases accurate predictions of airfoil 
pressure distributions can be obtained from an inviscid flow calculation 
without consideration of viscous boundary layer effects. If a boundary layer 
growth prediction is desired, the inviscid analysis can be combined in a 
non— interactive mode with a boundary layer analysis (e.g. Refs. 9 and 10) to 
predict both the pressure distribution and boundary layer development. The 
non-interactive inviscid flow-boundary layer calculation can give accurate 
predictions as long as viscous displacement effects are small and can even be 
used when limited regions of local separated flow are present (Ref. 11). 
However, for those cases in which viscous displacement effects significantly 
alter the inviscid pressure distribution, an alternate procedure is required. 

For the airfoil flow field problem significant viscous displacement 
effects are most pronounced in flows having regions of significant 
separation. In the presence of significant separation, the observed pressure 
distribution will differ considerably from that predicted from inviscid flow 
considerations. The actual pressure distribution corresponds to that around 
a body equivalent in shape to the airfoil plus a displacement correction (for 
viscous displacement effects), and in the presence of large separated regions 
the displacement correction is not small. In such cases an anlysis which is 
more complete than a purely inviscid analysis is required. One possibility 
for solving the separated airfoil flow field problem is the boundary layer 
strong interaction approach. In this approach an inviscid analysis and a 
boundary layer type analysis are solved so that the viscous displacement 
effects resulting from boundary layer growth influence the inviscid pressure 
distribution. Although this approach can give good results for some cases, 
it does have certain drawbacks. Usually, the approach requires an iteration 
between the two solutions and in the case of subsonic flow the iteration is a 
global one; e.g., the inviscid analysis is solved for a given displacement 
surface. The inviscid pressure distribution is then imposed upon the 
boundary layer equations and these equations are solved to predict the 
boundary layer development including a new displacement surface and the 
process is repeated until convergence occurs. This iteration process may be 
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difficult to converge under some circumstances, for example when large 
regions of separation occur or when the flow is transonic. Furthermore, 
assumptions may be required to treat the boundary layer equations in 
separated regions and normal pressure gradients must be assumed negligible in 
the viscous flow region. The approximations required in separated regions 
and the neglect of transverse pressure gradients in the regions where viscous 
effects are important may lead to serious inaccuracies for flow containing 
significant regions of separation. The drawbacks associated with boundary 
layer strong interaction techniques have led some investigators to seek an 
alternate means of predicting airfoil flow fields which solve the entire flow 
via a single set of viscous flow equations. 

Two general types of analyses which solve the viscous flow equations 
throughout the entire region of interest are in current use. These are the 
so-called T thin-shear layer 1 analysis and the Navier-Stokes analysis. The 
thin shear layer equations are an approximate form of the Navier-Stokes 
equations which contain all pressure and convective terms, but retain only 
some viscous terms. The viscous terms retained are those appropriate for a 
thin shear layer flow in which the shear layer is aligned with one of the 
computational coordinate directions. Use of the thin shear layer equations 
allows simultaneous calculation of pressure distribution and viscous and heat 
transfer effects without resorting to an interaction analysis and as such may 
be a valid approach for certain flow cases. However, since only part of the 
stress tensor is retained, the equations omit important terms for flows in 
which the streamlines vary significantly from computational coordinate lines 
such as in the case of flow separation. Although the thin shear layer 
equations are solved in less CPU time than are the Navier-Stokes equations, 
the saving is usually not more than twenty per cent. Therefore, thin shear 
layer equations do not have any major computer run time advantage over the 
full Navier-Stokes solution. Since one major item of interest in the present 
study focuses upon airfoils with significant regions of flow separation, only 
Navier-Stokes analyses shall be considered in the present discussions. 

The initial airfoil analyses based upon the Navier-Stokes equations 
considered incompressible laminar flow. Early examples are those of Mehta 
and Lavan (Ref. 12) and Lugt and Haussling (Ref. 13). Mehta and Lavan solved 
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a stream function vorticity formulation of the laminar incompressible 
Navier-Stokes equations to predict flow about an impulsively started 
airfoil and Lugt and Haussling utilized an incompressible stream 
function-vorticity approach to investigate flow about an abruptly started 
elliptical cylinder. More recent incompressible stream function-vorticity 
analyses have focused upon various aspects of the airfoil flow field 
problem. For example, Mehta (Ref. 14) used a numerical scheme considerably 
more efficient than that of Ref. 12 to solve incompressible laminar flow 
about an airfoil oscillating through incidence regimes in which stall 
occurs. Wu and Sampath (Ref. 15) and Wu, Sampath and Sankar (Ref. 16) 
applied the Wu-Thompson integro-differential formulation (Ref. 17) to both 
the impulsively started airfoil and the oscillating airfoil problem. In a 
similar vein Kinney and Cielak (Refs. 18 and 19) have investigated unsteady 
airfoil flow fields and Lugt and Haussling (Ref. 20) have investigated the 
time scale required to establish the Joukowski condition in incompressible 
flow. Finally, Thompson and his coworkers (e.g. Ref. 21) have calculated the 
flow about a variety of airfoil shapes and Hodge and Stone (Ref. 22) have 
investigated stalled airfoils using an incompressible primitive variable 
approach. 

Other investigations have considered compressible airfoil flow fields. 
Verhoff (Ref. 23) applied MacCormack's fully explicit method (Ref. 24) to the 
airfoil problem; however, since the procedure is fully explicit, a small time 
step is necessary to maintain numerical stability as a result of the locally 
refined mesh in the boundary layer and long computer run times result. In 
this regard conditionally stable schemes such as fully explicit schemes are 
not an optimum choice when mesh refinement is required for boundary layer 
definition; in these schemes the maximum allowable time step size is limited 
by the spatial step size leading to large run times. The time-step 
limitation problem, which is severe even in laminar flows, is magnified 
considerably in turbulent flows where a much finer spatial resolution is 
required in the boundary layer. On the other hand, unconditionally stable 
schemes (in a linear sense) such as some of the implicit schemes do not 
suffer from this characteristic. Both Deiwert's (Ref. 25) and Levy’s 
(Ref. 26) analyses are based upon MacCormack's hybrid implicit-explicit- 
characteristics scheme (Ref. 27). By virtue of an enlarged stability bound 
this new procedure is more efficient than the original MacCormack procedure 
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(Ref. 24) for airfoil calculations; however, it does present formidable 
coding problems. Implicit schemes, although more complicated to code than 
explicit schemes, do not present the formidable coding problems associated 
with the hybrid scheme. An implicit solution of the full Navier-Stokes 
equations has been developed by Gibeling, Shamroth and Eiseman (Ref. 28) who 
applied the Briley-McDonald (Ref. 29) numerical technique to the airfoil flow 
field. A similar approach has since been used by Sankar and Tassa (Ref. 30) 
to study an oscillating airfoil in a compressible low Reynolds number fluid. 
El Refaee , Wu and Likoudis (Ref. 31) have extended the integral approach of 
Wu and his associates (Refs. 15 and 16) to compressible flow. In another 
approach Steger (Ref. 32) used the thin shear layer equations in conjunction 
with the coordinate generation procedure of Thompson, Thames and Mastin 
(Ref. 33) to predict laminar flow about an airfoil. 

As expected most of the Navier- Stokes airfoil analyses were applied to 
laminar flow before turbulent flow calculations were attempted. In extending 
these analyses to the turbulent regime, several factors must be considered. 
Since turbulent flow requires considerably more resolution in the viscous 
wall boundary layer region than does laminar flow, a turbulent calculation 
requires very high near wall resolution and this high near wall resolution 
makes explicit procedures with their associated stability limits 
impractical. A second point concerns turbulence modelling. In general the 
airfoil flow field presents a very difficult turbulence modelling problem. 

The flow field contains regions of laminar, transitional and turbulent flow 
as well as significant separation regions. In cases where the flow is 
unsteady, shed vortices may be present. In general a turbulence energy model 
which includes a transition capability is required to analyze turbulent 
airfoil flow fields. However, useful information can also be gained by using 
a simpler model such as a mixing length model throughout. 

An early application to the airfoil turbulent flow problem was carried 
out by Shamroth and Gibeling (Ref. 34). The method was applied in Ref. 34 to 
airfoils at modest incidence. In addition it was applied to airfoils in 
stall by Shamroth and Gibeling (Ref. 35) and airfoils pitching at low 
incidence by Shamroth (Ref. 36). Other turbulent analyses are those of Tassa 
and Sankar (Ref. 37) and Sankar and Tang (Ref. 38) who studied a turbulent 
airfoil undergoing dynamic stall using an algebraic mixing length model. 
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The preseat report describes a work effort aimed at the turbulent 
airfoil flow field and contains results for a variety of steady and unsteady 
flow situations. Some of these have been previously reported in the open 
literature (Refs. 35 and 36), however, since the calculations presented in 
Refs. 35 and 36 were supported under the present contract, the results are 
repeated here. 
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ANALYSIS 

The Coordinate System 


The presence of bounding surfaces of a computational domain which do not 
fall upon coordinate lines presents significant difficulties for numerical 
techniques which solve the Navier-Stokes equations. If a bounding surface 
(such as the airfoil surface) does not coincide with a coordinate line, 
serious numerical errors may arise in the application of boundary conditions 
and considerable effort may be required to reduce these errors to an 
acceptable level. Although this problem arises in both viscous and inviscid 
flows, it is more severe in viscous flows where no-slip conditions on solid 
walls can combine with boundary condition truncation error to produce 
numerical solutions which are both qualitatively and quantitatively in 
error. Thus coordinate systems are sought in which each no-slip surface of 
the specific problem falls on a coordinate line. Such a system is termed a 
body-fitted coordinate system. Several approaches are available to form a 
body-fitted coordinate system. Among the coordinate system candidates are 
conformal coordinate systems such as that used by Mehta (Ref. 14), systems 
based upon solution of a Poisson equation such as those developed by Thompson 
and his coworkers (e.g.. Ref. 33) or Haussling (Ref. 39) and a constructive 
system. 

The approach used in the present effort is a constructive approach in 
which the required airfoil is by defintion a coordinate line and in which 
grid point placement is specified by the user. The procedure was developed 
originally for the isolated airfoil problem by Gibeling, Shamroth and Eiseraan 
(Ref. 28) and explained in general terms by Eiseman (Ref. 40). The 
coordinate system generated by the constructive process has several 
advantages. The system allows packing of grid points in regions where high 
grid resolution is required. In general, the high resolution regions are 
required near the airfoil surface (where the boundary layer is found) and in 
the vicinity of the airfoil leading edge where rapid streamwise changes are 
present. In addition, although the grid has a branch cut emanating from the 
airfoil trailing edge, metric data are continuous across the branch cut. 
Furthermore, although the grid is nonorthogonal, the amount of 
nonorthogonality is not large. Finally, as applied to the airfoil problem 
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the metric data remain smooth from grid point to grid point. A computer 
generated plot of the airfoil coordinate system is given in Fig. 1 where for 
clarity not all lines are shown. The actual grid used has very high 
resolution near the surface where the transverse grid spacing is of the order 
of 10 -5 chords and near the leading edge where the streamwise grid spacing 
is of the order of 10“3 chords. 

The coordinate system consists of a set of two families of curves; the 
£ = constant curves such as line HI in Fig. 1 and the n = constant curves 
such as ABCD or A'ED' in Fig. 1. The coordinate system is constructed by 
first forming the inner loop A’ED' which includes the airfoil. The airfoil 
may be either specified by an analytic equation or by discrete data points. 

If an equation is used, then construction of the inner loop is 
straight-forward. If the airfoil is specified by discrete data points, then, 
in general, the points required on the inner loop will not coincide with any 
point used for airfoil specification. In this case, a curve fit is used to 
obtain the required inner loop points. The curve fit used is based upon a 
local parabolic fit. For any given point required on the inner loop, a 
parabola is fit through three adjacent specifying points, two on the right 
and one on the left. A second parabola is then fit through the two points on 
the left and one on the right. The location of the required point is 
obtained via a weighted average of these curve fits with the weighting factor 
being determined by the distance from the required point to the center 
specifying point of each parabola. Alternatively the airfoil may be fit in a 
piecewise manner through a least squares polynominal constrained to maintain 
continuous first and second derivatives at joining points. 

This is followed by constructing an outer loop ABCD which consists of 
line segments AB and CD at a specified distance from the airfoil chord line 
and a frontal curve BC. Both the inner and outer loops are then represented 
by parametric curves 


x = x (s) , y = y(s) (1) 

where the parameter varies from zero to unity. The present coordinate 
generation process utilizes a multi— part transformation. First x and y are 
expressed as a function of s', the physical distance along the curve. Then 
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s f is related to s via a hyperbolic tangent parameterization centered about 
the leading edge for the inner loop and a cubic polynomial representation for 
the outer loop. The inner loop transformation parametric representation is 
chosen so as to have rapid variation in the airfoil leading edge region. 

Both the inner loop hyperbolic tangent transformation and the outer loop 
cubic transformation are applied between s '/2 and 0 and between s ’/2 and s’ 
on each loop. This ensures that corresponding points on the branch cut will 
be equi-distant in s from the branch cut end points, points A T and D T . This 
property is required if corresponding branch cut points are to fall on the 
same pseudo-radial line. 

Having specified loop 1, the inner loop, and loop 4, the outer loop, the 
construction process now specifies loop 2 and loop 3. Both loop 2 and loop 3 
are located between loop 1 and loop 4 with loop 2 being a normal distance A 2 
from loop 1 and with loop 3 being a normal distance A 3 from loop 4. Loops 2 

It M 

and 3 have parameters S 2 and S 3 associated with them and the method of 

determination of these parameters is discussed subsequently. Before 

discussing this point it is necessary to introduce a pseudo-radial parameter 

r and a position vector P associated with each loop. The radial parameter is 

defined at the downstream boundary (line A'A) as the distance from the loop 

in question to the inner loop divided by the distance from the outer loop to 

the inner loop. Thus ri =0, r 2 = A 2/ R MAX> r 3 = ( R MAX “ a 3>/RmAX and 

r 4 = 1 where R^AX is the distance from the inner loop to the outer loop. 

*► 

The position vector P-^(t) is a vector whose components are the x and y 
coordinates of the i^h loop (i = 1 , 2, 3, 4) when the parameter = t 
where t is a number between 0 and 1. With the definition of these quantities 
it is possible to introduce the general position vector P(r,t) where 

P(r,t) = ( I - r ) 2 ( I - a ( r)P ( (t) + ( + 2)( i - r) 2 rP g ( t ) 

/ n \ 

+ r 2 [l-a 2 (i-r)]p 4 (t) + (a 2 + 2)r 2 (l-r)P 3 (t) 


a 


i 



2 

° 2 = 3(1 - r 2 )- | 


(3) 


10 



It should be noted that at r = 0, P(0,t) = Pi(t) and at r = 1 , 
P(l,t) = P4 ( t ) • Further since at r = 0, 


dP 

17 


(0,t) = [p 2 (t) - t )](a, + 2) 


( 4 ) 


and at r = 1 


dp 

dr 


(l,t) = [P 4 (t) - P 3 (t)](a 2 + 2) 


(5) 


specification of the derivatives at the inner and outer boundaries determines 
the parametric representation of intermediate loops 2 and 3. Since 3?/3r 
implies 3x/3r and 3y/3r, specification of 3P/ 3r at the boundaries controls 
the angle which the pseudo-radial coordinate line forms with the boundary 
line. Thus, the four loop method allows specification of the boundary point 
locations and coordinate angles at these boundaries. 

After loops 2 and 3 are parameterized to satisfy the coordinate angle at 
the boundary points, the grid is constructed as follows. If the grid is to 
contain M pseudo— radial lines (such as line HI of Fig. 1) and N 
psuedo-azimuthal lines (such as line QPR) , the values of the pseudo-radial 
coordinate are 


r(i) = i /( N - 1 ) i= 0,1 ,2, .... N-l 


and the values of the pseudo-azimuthal coordinate are 


t(j) = j /( M-l) j = 0,1,2, ... , M-l 


Then the position vector for each point in the grid is given by Eq. (1). 

The preceding has assumed a uniform spacing in the radial direction. If 
radial grid point concentration is desired, it is simply necessary to assume 
a radial distribution function. The present analysis assumed a distribution 
function 


tanhD(l-r) 

tanhD 


( 6 ) 
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which concentrates points in the wall region. Grid points then were chosen 
at r(i) = (i)/(N-l) and the analysis proceeded as outlined. Once the grid 
point locations are obtained, the required metric data can be calculated by 
numerical differentiation. 


Mean Flow Equations 

A solution of the compressible, time-dependent Navier-Stokes equations 
in conjunction with a suitable turbulence model would serve to predict the 
flow field for both laminar and turbulent flows. The form of the equations 
expressed in the more common coordinate systems can be found in standard 
fluid dynamic texts and the equations themselves have been derived in general 
tensor form by McVittie (Ref. 41) for inviscid flow and by Walkden (Ref. 42) 
for viscous flow. 

One possible approach for solving the equations in general nonorthogonal 
form is the strong conservation approach such as that used by Thomas and 
Lombard (Ref. 43). A second possible approach solves a set of equations in 
which the metric coefficients do not appear within derivatives (quasilinear 
form). In both cases the independent spatial variables are transformed from 
the Cartesian coordinates (x,y) to a new set of coordinates (£,q) where 


£ = £ (x,y,t) 

V = 

T-t (7) 


The strong 


conservation form of the equations then becomes 
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where 


0 ' My - Mx 


w * 



, F 



/JU 2 + P ) , 

P UV / 


l ^ \ 

G ' (/“Ip ) 




( 9 ) 


The quasilinear form of the equations is expressed as 


dW aw 

TT + e t ae" 

i 

Re 


dF dG 

+ f «dT + f ylT 

dF, aF, 

**nr + 


aw aF dG 

+ drj Qjj * 


+ € 


dGj 

ylf 


dG, ' 
)f~drf 


(10) 


The problem of proper equation form in non-Cartesian spatial variables 
has been discussed by several investigators (e.g., Refs. 34, 43 and 44). If 
the strong conservation form of the equations is to be used, then care must 
be taken to evaluate the metric data by a method which is consistent with a 
control volume approach (Ref. 43). Usually this requires numerical 
evaluation of the metric data even if an analytic functional relationship for 
the transformation is available. The analytic representation of the metric 
data, 5 X > 5y, etc., when combined with the strong conservation form of 
the equations leads to significant error for as straightforward a calculation 
as low Reynolds number flow about a circular cylinder (Ref. 34). The present 
effort utilizes the quasilinear form of the equations since this form is much 
less sensitive to the form of metric evaluation and gives good results for 
both numerical and analytic evaluations of the metric data. Furthermore, as 
shown in Refs. 45-47, the quasilinear equations have been used with good 
results for transonic flow in a cascade environment. 

A final item related to the choice of equations is the choice of 
dependent variables. In the present approach the density and velocity 
components are used as dependent variables. The energy equation is replaced 
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by an assumption of constant total temperature thus leading to a relation 
between pressure, density and velocity. 

u 2 + w 2 

P = — — ) (ID 

where R is the gas constant, T° is total temperature and C p is specific 
heat. 

It should be noted that the energy equation can be solved with the 
momenta and continuity equations at the cost of adding an additional 
governing equation which increases computer run time. Calculations of this 
type in transonic cascades which include comparison with heat transfer data 
have been made by Weinberg, Yang, Shamroth and McDonald (Ref. 48). For 
steady airfoil flow fields this assumption is reasonable. For unsteady flow, 
it represents an approximation as can be noted from examination of the 
unsteady total temperature equation. However, as discussed in Ref. 49, for 
the cases considered here this assumption should still be valid. 

The Turbulence Model 

Since the present effort is concerned with high Reynolds number 
turbulent flows, it is necessary to specify a turbulence model. The results 
presented were obtained primarily with a mixing length model. The mixing 
length model assumes the existence of a mixing length, Jt, and then relates an 
eddy viscosity, Uj, to the mixing length by 


P-T 



( 12 ) 


For flow regions upstream of the leading edge where the flow is attached the 
mixing length is determined by the usual boundary layer formulation 

i=*yo 03) 
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where < is the von-Karman constant, D is a sublayer damping factor and £ raax 
is taken 0.09 6 where S is the boundary layer thickness. The damping factor, 
D, which has for the most part been utilized is the van Driest damping factor 

-V + /27 

D = (I - e Y ) (14) 

where y + is the dimensionless coordinate normal to the wall, yu T / v# 

When the mixing length formulation is used in a boundary layer 
environment, 6 is usually taken as the location where u/u e = 0.99. 

However, this definition assumes the existence of an outer portion of the 
flow where u e is independent of distance from the wall and assumes that the 
location where u e becomes independent of distance from the wall marks the 
end of the viscous region. In an airfoil Navier— Stokes calculation no such 
clear flow division occurs as u approaches the upstream velocity, u», as 
distance from the wall increases. Therefore, the boundary layer thickness, 
5, is set by first determining u raax , the maximum velocity at each given 
streamwise station and then setting 6 by 

8 ‘ 2 - 0 y t“/w k .) <15) 

i.e., 6 is taken as twice the distance from the wall to the location where 
u/u m ax = ki . Two values of kj were used; these were 0.80 and 0.90. If the 
flow is separated, then a minimum mixing length is set by 

L in = 0.1 hD (16) 

mm 

where h is the local height of the separated region. In the wake the mixing 

length is made proportional to the wake thickness, 5, and a linear growth of 

6 with distance is assumed based upon classical free jet boundary growth 
(e.g.. Ref. 50). With this assumption 

S = (S ps + S ss )+ 0.2 (X - X TE ) (17) 

where 6pg and SgS are the pressure and suction surface trailing edge 
boundary layer thicknesses and is the trailing edge location. 
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The mixing length, A, is taken as 0*26, The viscosity is smoothed between 
regions obtained using the wall formulation for H and the wake formulation 
for H. Having obtained the turbulent viscosity, y T , the turbulent stress, 
-pui'uj' is given by 





au^ 

a *k 



(18) 


Although the mixing length model does not include a transition model, 
transition can be simulated by specifying a location upstream of which the 
flow is laminar. This corresponds to forced transition. Even if no forced 
transition is assumed, the flow in the leading edge region will be laminar as 
the boundary layer thickness becomes very small in this region. 

A second turbulence model which was implemented for a case of an 
NACA 0012 airfoil at 6° incidence was the two-equation k-e model. This model 
is well known (e.g. Refs. 51—54), and has been used by several investigators 
(e.g. , Ref. 48). In brief, the model is based upon a turbulence energy 
equation 


dpk d/auk 

at + ax 


+ 



auj_ 

ax k 


+ 



(19) 


and a turbulence dissipation equation 


dpt dpue dpvc 

~dT + + ~dy~ “ 

« / au, au„ \ du-. 

+ C| k dx k + dx. j ax k 



( 20 ) 


The turbulence viscosity is then obtained via the Prandtl-Kolmogorov relation 

H-j ~ />C^k 2 /< f(y/8) (21) 
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where Cp is a turbulence structural coefficient and f(y/ <$) is a factor used 
to ensure small turbulent viscosities at locations far from the airfoil. 

The function f(y/<$) is taken as 

f(y/8) =1.0 y £ S 

f(y/8) rf-^; 1 ' 01 y > 8 (22) 


where b is a constant. 

In the present analysis the following values were assumed 

a~ € “ 13 
o- k - 1.0 
C, = 1.55 

and C y and C 2 were made functions of turbulence Reynolds number, 
R x = pk 2 /ue. 


(23) 


* 0-09 exp [-2 5(1. + R r /50.)] (24) 

C 2 - 2.o{l O-0 3[ exp(-R x )]} (25) 


Numerical Procedure 

The numerical procedure used to solve the governing equations is a 
consistently split linearized block implicit (LBI) scheme originally 
developed by Briley and McDonald (Ref. 29). A conceptually similar scheme 
has been developed for two-dimensional MHD problems by Lindemuth and Killeen 
(Ref. 55). The procedure is discussed in detail in Refs. 29 and 56. The 
method can be briefly outlined as follows: the governing equations are 

replaced by an implicit time difference approximation, optionally a backward 
difference or Crank-Nicolson scheme. Terms involving nonlinearities at the 
implicit time level are linearized by Taylor expansion in time about the 
solution at the known time level, and spatial difference approximations are 
introduced. The result is a system of multidimensional coupled (but linear) 
difference equations for the dependent variables at the unknown or implicit 
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time level. To solve these difference equations, the Douglas-Gunn (Ref. 57) 
procedure for generating alternating direction implicit (ADI) schemes as 
perturbations of fundamental implicit difference schemes is introduced in its 
natural extension to systems of partial differential equations. This 
technique leads to systems of coupled linear difference equations having 
narrow block-banded matrix structures which can be solved efficiently by 
standard block-elimination methods. 

The method centers around the use of a formal linearization technique 
adapted for the integration of initial-value problems. The linearization 
technique, which requires an implicit solution procedure, permits the 
solution of coupled nonlinear equations in one space dimension (to the 
requisite degree of accuracy) by a one-step noniterative scheme. Since no 
iteration is required to compute the solution for a single time step, and 
since only moderate effort is required for solution of the implicit 
difference equations, the method is computationally efficient; this 
efficiency is retained for multidimensional problems by using what might be 
termed block ADI techniques. The method is also economical in terms of 
computer storage, and in its present form requires only two time levels of 
storage for each dependent variable. Furthermore, the block ADI technique 
reduces multidimensional problems to sequences of calculations which are 
one-dimensional in the sense that easily solved narrow block-banded matrices 
associated with one-dimensional rows of grid points are produced. A more 
detailed discussion of the solution procedure is discussed by Briley, Buggeln 
and McDonald (Ref. 58) and is given in the Appendix. 


Boundary Conditions 

An important component of the airfoil analysis concerns specification of 
boundary conditions. The present analysis requires boundary conditions to be 
set along the lines, £ = 5min> £ = ?max» ^ = Rmin ar *d B With 

the coordinate system sketched in Fig. 1, £ = (line AA f ) and 

5 = 5max (line DD T ) are downstream boundaries. In the early work done 
under this effort derivatives were set to zero at this boundary and function 
conditions specified on the remainder of the outer boundary. On the airfoil 
surface no-slip conditions are used in conjunction with an inviscid momentum 
equation (which for no motion and no heat transfer reduced to zero density 
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gradient) as boundary conditions and either the turbulence energy or its 
derivative was specified at the surface when turbulence energy was included 
as a variable. More recently the boundary conditions were modified based 
upon a suggestion by Briley and McDonald (Ref. 59). Following this 
suggestion, static pressure is specified along with velocity derivatives 
along the downstream boundaries (lines AA’ and DD’) and along the aft portion 
of the outer boundary (line segments AB and CD). Total pressure, angle of 
incidence and the density derivative are specified along the outer boundary 
segment BC. This approach was used successfully (Refs. 45-48, 60) in a 
Navier-Stokes solution to the cascade problem and has been incorporated into 
the airfoil analysis. 

In addition, calculations have been made with the full transverse 
momentum equation rather than the normal pressure gradient equal to zero as a 
wall boundary condition. Little difference was noted in the solution 
although the full transverse momentum equation boundary condition 
occasionally showed some tendency toward numerical instability and may 
require smaller time-steps. Finally, calculations have been made in which 
tunnel wall boundary conditions are simulated by specifying the flow 
direction and a full slip condition along AB and CD. 

Artificial Dissipation 


One major problem to be overcome in calculating high Reynolds number 
flows using the Navier-Stokes equations is the appearance of spatial 
oscillations associated with the so-called central difference problem. When 
spatial derivatives are represented by central differences, high Reynolds 
number flows can exhibit a saw tooth type oscillation unless some mechanism 
is added to the equations to suppress their appearance. This dissipation 
mechanism can be added implicitly to the equations via the spatial difference 
molecule (e.g., one-sided differencing) or explicitly through addition of a 
specific term. The present author favors this latter approach for two 
reasons. First, if a specific artificial dissipation term is added to the 
equations, it is clear precisely what approximation is being made. Secondly, 
if a specific term is added to suppress oscillations, the amount of 
artificial dissipation added to the equations can be easily controlled in 
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magnitude and location so as to add the minimum amount necessary to suppress 
spatial oscillations. 

The results presented herein basically utilize two levels of artificial 

dissipation. During initial phases of the present study a term of the form 
2 2 

v art ^ <J>/9z was added to the governing equation where <j> = p, u, v for the 
continuity, x-moraentum and y-momentum equations respectively and v art is 
determined by 

U z AZ I 

v + (l/ art^ z £ CT Z < 26 > 

In the above equation AZ is the distance between grid points in a given 
coordinate direction, Ug is the velocity in this direction, 0 % is the 
artificial dissipation parameter for this direction and v is the kinematic 
viscosity. The equation determines v art with v art taken as the smallest 
non-negative value which will satisfy the expression. It should be noted 
that in two space dimensions each equation contains two artificial 
dissipation terms, one in each coordinate direction. For example, the 
streamwise momentum equation expressed in Cartesian coordinates would contain 
the artificial dissipation terms 

/ \ . / \ ^ 2 U 

-^5 (27) 

The parameter o z was taken as 0.5. 

During the time period of the present contract effort, various methods 
of adding artificial dissipation were investigated in Ref. 46 and these were 
evaluated in the context of a one-dimensional model problem. The model 
problem used was one-dimensional flow with heat transfer. Flow was subsonic 
at the upstream boundary, accelerated via heat sources until a Mach number of 
unity was reached and then accelerated by heat sinks. The exit back pressure 
was raised to cause a shock to appear in the supersonic region. This basic 
one— dimensional problem contained many relevant features including strong 
accelerations and, therefore, it served as a good test case for evaluating 
various forms of artificial dissipation which could be used. Several 
different types of artificial dissipation terms were considered, and 
it was concluded that for the present numerical method, a second order 
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artificial dissipation approach such as represented by Eqs. (26) and (27) is 
suitable. However, instead of the parameter, o, being 0.5 it should be set 
at approximately 0.05. At this value, sufficient artificial dissipation is 
added to the equations to suppress spurious oscillations, but the amount 
added does not significantly change the physical solution. Furthermore, as 
shown in Refs. 47-48, the solutions obtained with a = 0.050 change only 
little when a is lowered to 0.025. Other confirmations of this approach were 
given in Refs. 61 and 62 as well as in this present effort. Based upon these 
studies, it has been concluded that a < 0.1 gives accurate representation for 
most two-dimensional flows. 

An alternate method of including artificial viscosity is the so-called 
conservation form in which dissipation terras of the type 

sr {(■'■"•h f^} + 77 t (■'"«)» 77} < 28 > 

are added to the equations. This form confines the integrated effect of the 
terms over the computational domain to the computational boundaries since 

ffix{ ( l '« r *)x fr }^ d y- 

A 

where A is the area of the computational domain and S is the boundary line. 

If the artificial viscosity and its derivative are set to zero on the 
boundary, global conservation is obtained. However, local conservation 
obviously is not. Calculations have been run with this latter form of 
artificial dissipation and for the cases run no significant differences 
appear in the predicted flow fields due to artificial dissipation of the form 
given in (28) as opposed to that given in equation (27). 


Convergence and Run Times 

When a steady flow is sought via a time marching technique, the question 
arises as to when convergence is obtained. In considering this question, 
several factors must be taken Into account* First of all, not all flows 
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reach a steady state. For example, airfoils at high incidence which shed 
vortices or airfoil flows in which a shock wave is present may never become 
truly steady. In the former case, vortices are shed in some quasi-periodic 
manner and the unsteadiness has a large time scale. In the latter case, the 
shock position may move leading to an unsteadiness with a small time scale. 
Obviously, in these cases no steady flow solution is guaranteed. Secondly, 
the numerical technique used may hinder complete convergence. For example, 
in the present approach the turbulent viscosity is lagged by one step in time 
and this interaction between the viscosity evaluation and the mean flow 
calculation may hinder or even prevent complete convergence. 

During the early part of the present effort when steady state solutions 
were sought, the calculated results were monitored to assess convergence. In 
particular, the surface pressure distribution, the location of separation 
points and the velocity field in the vicinity of the airfoil were monitored 
and when they ceased to vary significantly for significant time scales, the 
calculation was terminated. During latter parts of the effort, an additional 
criterion was added. This was based upon the evaluation of the residual for 
each equation. The residual of each equation is obtained by setting the 
time-derivative term to zero, and placing all remaining terms on the 
right-hand side of the equation. The sum of all terras on the right hand side 
of the equation defines the residual. Obviously, when the residual is zero 
the equations satisfy a steady state solution. Both the maximum residual 
throughout the domain for each equation and the average residual within the 
domain for each equation were monitored. These usually could be decreased by 
between two to four orders of magnitude during the run when steady solutions 
were sought. 

However, even the presence of residuals requires interpretation. As 
previously discussed, these could be indicative of flow unsteadiness. Also, 
relatively large residuals occur at the airfoil cusp trailing edge. In 
general, calculations for which steady solutions are sought are initiated 
from a very simple initial flow field. The initial flow field has constant 
pressure throughout and a velocity field identical to that at upstream 
infinity with a simple boundary layer correction. For steady flows converged 
solutions are usually obtained within 70 time steps. 

In regard to run time, the current code is a general research type code 
which was created with flexibility in mind. Therefore, items such as the 
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equations solved, boundary conditions specified, dependent variable choice, 
form of artificial dissipation, etc. can be changed with only a small amount 
of effort. Obviously, the price of flexibility is increased computer run 
time. In the present code the run time for a 141 x 39 grid is approximately 
15 cpu secs per time step on a CYBER 203. The code used is not fully 
optimized for scalar operation and has no vectorization. It is estimated 
that the run time could be reduced to 7 sec per time step with further scalar 
optimization without compromising existing generality and an additional 
factor of ten could be obtained through vectorization. 
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RESULTS 


Results obtained under this contract effort include both steady and 
unsteady airfoil calculations under a variety of flow conditions. Early 
results appear in Ref. 34. Results obtained since Ref. 34 are presented in 
the present section. 

Low Mach Number - High Incidence NACA 0012 Airfoil 

The low Mach number - high incidence calculation was made for an 
NACA 0012 airfoil immersed in a free stream having a Reynolds number of 10 6 
and a Mach number of 0.148. The calculation was run prior to the artificial 
dissipation study (Refs. 47 and 48) and consequently a dissipation parameter, 
o, equal to 0.5 was used. Since the flow was subsonic, the results are 
expected to be qualitatively correct with more numerical dissipation than 
actually required; the major discrepancy is expected to be in the leading 
edge suction peak region. Results of this case which have been presented in 
E- e f» 35 are reproduced in a condensed form here. 

The calculation was run on a 'C' type grid having 81 pseudo-radial lines 
and 39 pseudo-azimuthal lines. The grid was highly non-uniform with the 
first point off the airfoil a distance of 0.2 x lO - ** chords away from the 
airfoil. In contrast, the last point was placed four chords away from the 
airfoil with a radial spacing of 0.5 chords. Similarly, the streamwise 
(pseudo-azimuthal) grid was concentrated in the vicinity of the airfoil 
leading edge with the minimum spacing being approximately 0.4 x 10~ 2 chords. 
The calculation was initiated from a converged solution for the airfoil at 
6 degrees incidence. The incidence was then changed to 19 degrees via the 
equation 


a 


a o + 



COS U) (t - t 0 


>] 


t Q < f < (t 0 + 7r ) /CO (30) 
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where 


a 0 =6°, Aa = 13°, cu = 5 , t Q = 1.20 (31) 

For t>(t 0 +II)/a>, the incidence was held constant at 

a = a 0 + Aa (32) 

It should be noted that the dimensionless frequency, k = toc/2Ueo = 2.5, 
represents a very high value for airfoil calculations with the 13° ramp 
amplitude. 

The results of the calculation during the ramping period are presented 
in Figs. 2 and 3. Figure 2 shows the pressure coefficient distribution at 
various incidence angles. At six degrees the pressure distribution is 
typical of that found for a steady airfoil; the suction peak has been smeared 
and diminished due to insufficient streamwise resolution and the relatively 
large amount of artificial dissipation as discussed previously. As the 
incidence changes from 6 to 9 degrees the rapid motion, particularly in the 
trailing edge region, causes high pressure to appear on the lower side of the 
airfoil and low pressures to appear on the upper side. It should be noted 
that the velocity of the airfoil trailing edge relative to the inertial frame 
reaches a maximum value of 0.4 U<» and, therefore, large deviations from the 
steady solution are to be expected. The situation becomes more pronounced at 
12.5 degrees; however, by 14 degrees a tendency to return to the usual static 
airfoil pressure distribution appears. Finally, at the last incidence angle, 
19 degrees, (t = 1.93), the basic pressure distribution is approaching the 
type expected for a steady airfoil with no evidence of stall. The location 
of the separation points is presented in Fig. 3. At the initiation of the 
ramp motion no separated flow was present; however, separation appeared soon 
after the ramp motion began and the trailing edge separation point moves 
continuously upstream as shown in the figure. During this process the 
separated region remains very thin and has only a minimum viscous 
displacement effect upon the outer nominally inviscid flow. 

After cessation of the motion, the flow continues to develop and the 
pressure distribution undergoes radical changes as shown in Figs. 4-6. The 
major changes occur in the airfoil leading edge region where the suction peak 
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appearing on the airfoil upper surface continues to drop in magnitude from a 
value of approximately 6.8 at t = 1.83 (just after the cessation of airfoil 
motion) to a value of approximately 1.2 at t = 5.38. A unit increment in t 
represents the time required for a particle moving at free stream velocity to 
transverse a distance of one chord. The drop in the suction peak and the 
accompanying decrease in airfoil lift exhibited in Figs. 4-6 are consistent 
with the development of airfoil stall. The calculation also predicts a minor 
movement of the airfoil front stagnation point towards the geometric leading 
edge. In addition to the loss of lift, the analysis predicts a pressure 
perturbation to initiate at t * 3.7 (see Fig, 5) and then move downstream at 
a speed of approximately 35% free stream velocity. Although quantitative 
comparisons between this prediction and data are not available, the predicted 
flow seems physically realistic. 

Upon reaching 19 degrees, the motion ceased and the airfoil flow field 
was allowed to develop at 19 °. A comparison of the calculated results and 
the measured data of Young, Meyers and Hoad (Ref. 63) for an airfoil at 19.4° 
incidence is presented in Fig. 7, Figure 7 compares the predicted and 
measured values of the zero velocity line. Below this line, the flow is 
directed toward the leading edge and above this line the flow is directed 
toward the trailing edge. The predicted values are shown as a function of 
time. During the ramping process, the separated region present was too thin 
to be shown on the scale of Fig. 7 and the results shown are at times well 
past the cessation of the ramping motion at t a 1.9. The results presented 
in Fig. 7 show the growth of the backflow velocity zone with time, and at the 
latter times shown the backflow zone position has converged over most of the 
airofil as continued growth is confined to regions in the vicinity of the 
airfoil trailing edge. As can be seen, the comparison between the calculated 
zone location and that measured by Young, Meyers and Hoad is very 
good. 

A vector plot of the velocity field as measured by Young, Meyers and 
Hoad is shown in Fig. 8. These results show a large separated region to be 
present over the airfoil upper surface with separation initiating in the 
immediate vicinity of the airfoil leading edge. A vortex appears to be 
centered at roughly the eighty percent chord location. The data (not shown 
on this figure) indicated that the wake closure point was located well 
downstream of the airfoil trailing edge and above the airfoil suction 
surface. 
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Another feature is the appearance of a very strong shear layer in the airfoil 
trailing edge vicinity where the suction surface and pressure surface flow 
fields meet. Finally, the calculated results. Figs. 9 and 10, indicate that 
flow is entrained into the recirculation region from two sources. One source 
is the flow region above the recirculation zone. The second source is the 
flow which originates on the airfoil pressure surface, then passes into the 
mixing layer which forms at the airfoil trailing edge and finally is 
entrained into the recirculation region from below. 

Predicted velocity vector fields are shown in Figs. 9 and 10. These 
figures represent the flow field at times t x and ti+At where At is the time 
required for a free stream particle to move a distance of one chord length. 

The analysis predicts the formation of a large separation region which 
initiates very near the airfoil leading edge; the calculated flow field is in 
qualitative agreement with the data shown in Fig. 8. Other similarities 
between data and calculation can be found in the vortex formation and in the 
strong shear layer which appears at the airfoil trailing edge. In addition, 
the calculated flow field was characterized by significant flow unsteadiness 
in the leading edge region which limited the permissible maximum time step. 
This characteristic of unsteady leading edge flow also appeared in the 
experimental s tudy • 

In regard to other features, the analysis showed the vortex to be moving 
downsteam at a velocity of approximately 0.3 Ua>; however, no regular 
shedding pattern was observed in the experiment. Some comments on this are 
in order. First of all, the calculation was not run long enough to determine 
if a regular shedding will result although the first vortex being formed 
definitely appears to be in the process of shedding. Secondly, although the 
experiment did not detect any regular shedding pattern, it is possible that 
an irregular shedding did occur. 

Calculated vorticity contours at the two times are shown in Figs. 11 
and 12. The vorticity contours presented correspond to normalized values of 
-100, -25,-10, -5, 0, 5, 10, 25, 100. In both figures, the vorticity on the 
airfoil pressure surface is confined to the boundary layer whereas that on 
the suction surfaces occurs in two locations. One region of vorticity is 
located in the wall layer close to the airfoil surface; the second region is 
a 'tongue-like* region extending from the vicinity of the airfoil leading 
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edge into the 'free stream'. This contour line represented by the value 5 is 
a region of a local maximum vorticity. As can be seen by comparing Figs. 11 
and 12, the tongue-like region of vorticity appears to break off and be 
convected downstream as a local concentrated region (See Fig. 12). This may 
be interpreted as the initiation of a shed vortex. A third area of high 
vorticity concentration occurs at the airfoil trailing edge where the sharp 
mixing layer is present. 

A closer examination of the predicted flow field shows the emergence of 
an inner counter-clockwise rotating separation zone which occurs under the 
main suction surface separation zone. As can be seen in Figs. 8-10, the 
major separated region is a large region of clockwise rotation. However, a 
detailed vector plot of the mid-chord portion of the suction surface 
presented in Fig. 13 shows a secondary separation region of counter-clockwise 
rotation completely embedded within the primary separated zone. The 
stagnation point location, the flow separation at the stagnation point, the 
acceleration about the leading edge and the initiation of flow separation are 
all shown clearly. A detail of the leading edge region is shown in Fig. 14. 
Static pressure contours are presented in Ref. 35. 
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0.5 Mach Number - High Incidence NACA 0012 Airfoil 


A calculation similar to the low Mach number - high incidence case 
discussed previously was made for an NACA 0012 airfoil immersed in a free 

g 

stream of Reynolds number 1.4 x 10 and Mach number of 0.5. The calculation 
was initiated from 6 degrees incidence with a cosine ramp motion, 

Eqs. (30-32), bringing the airfoil to 19 degrees incidence. In general, the 
results for the 19°, Moo = .5 case were similar to the results obtained for 
the 19°, - • 15 case discussed previously. The flow did not become 

steady during the calculation as a large separation zone appeared over the 
suction surface. A velocity vector plot showing the flow field at a time 
2 T re f past the cessation of airfoil motion is shown in Fig. 15; T re f is 
the time for a free stream particle to move a distance of one chord. Several 
features of the flow are clearly evident in this figure. The flow approaches 
the airfoil leading edge and branches with part of it passing over the 
suction surface and part passing over the pressure surface. A detail of the 
leading edge region is shown in Fig. 16. Figures 15 and 16 clearly show the 
flow branching and the acceleration about the airfoil leading edge. The flow 
on the pressure surface remains well-behaved with a thin boundary layer 
developing on the airfoil surface. The situation on the suction surface is 
much different with separation occurring almost at the airfoil leading edge 
and a large separation zone being formed over the suction surface. The 
separation zone and accompanying vortex is also shown clearly in Fig. 17 
which gives vorticity contours, the vorticity being defined as V x V. The 

vorticity pattern is generally similar to that obtained for the Mo = .15 
calculation with the maximum core value in the present case being 
approximately twice that of the previous case* Static pressure contours are 
presented in Fig. 18. 

Flow About an Airfoil Oscillating in Pitch at Low Incidence 

The first oscillating airfoil case considered was flow about an 
NACA 0012 airfoil immersed in a free stream at a Reynolds number based upon 
airfoil chord of 0.26 x 10 7 , and a Mach number of 0.20. The case was run as 
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initial demonstration case for an oscillating ai oil prior to attacking 
the dynamic stall problem. The airfoil was assumed o oscillate in pitch 
about its quarter chord point with the incidence giv *n by 


Ad r 

a = a Q + — [ 1.0 - cos u>(t - t 0 )J 


The reduced frequency, k = wc/2Uoo, was taken to be 0.12. The calculation 
was run with the dissipation parameter, o, set equal to 0.5 and was initiated 
from a steady solution at zero degree incidence and then followed the 
incidence variation given by Eq. (33). 

The predicted lift versus incidence curve is presented along with the 
steady and unsteady data of Grey and Liiva (Ref. 64) in Fig. 19. Considering 
first the experimental data, the unsteady curve shows a hysteresis loop. 
Furthermore, the general slope of the curve is less than that of the steady 
data and the unsteady lift at zero incidence is higher than that of the 
steady data (which is zero). The prediction shows the same general 
characterstics. The calculation was initiated at zero degrees incidence from 
a steady calculation and followed the theoretical quasi-steady lift-incidence 
curve until a « 4°. After reaching 4°, the lift predicted is less than the 
inviscid value and this is primarily a result of the under prediction of the 
suction peak resulting from the specification of a = 0.5 where a is the 
artificial dissipation parameter. The dissipation factor was taken as 
a = 0.5 since this calculation was made prior to the dissipation parameter 
study of Refs. 45 and 46. Upon reaching the maximum incidence, a = 10.5°, 
the curve forms a hysteresis loop as incidence decreases. This loop is 
somewhat more pronounced than that measured. After reaching the minimum 
value of a = 0°, the lift increases with incidence and at the last time 
calculated the loop is closing. Although the thickness of the predicted 
hysteresis loop is somewhat greater than that of the measured loop, the 
average slopes agree. In addition, both prediction and data show significant 
lift at zero incidence; this is in contrast to the quasi-steady calculation. 

A detailed examination of the flow field prediction shows the major 
contribution to the finite thickness lift loop results from the suction 
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surface boundary layer thickness for o < 0 being greater than that for a > 0 
at the same value of a. This result represents a lag in the boundary layer 
reaction of the pressure gradient which modifies the mid-chord and trailing 
pressure distribution. The mid-chord and trailing edge effect is somewhat 
modified by differences in the leading edge where the suction peak for a < 0 
is more pronounced than that for ot > 0. It should be noted that the loop 
calculation is a very sensitive one and its formation results from relatively 
small pressure changes on both the pressure and suction surfaces. Further, 
this interpretation should be regarded as approximate due to the value of 
artificial dissipation used. Velocity vector plots are given in 
Figs. 20-22. These figures clearly show the general flow pattern which 
includes the approach to the leading edge stagnation point, acceleration 
ground the leading edge and the boundary layer and wake development. 

A comparison of Figs. 20 and 22 shows that during the upstroke (ot > 0) the 
flow along the aft portion of the airfoil tends to align with the pressure 
(lower) surface. Furthermore, the differences in the suction surface 
boundary layer thickness and wake position are evident. Vector plots in the 
vicinity of the leading edge are shown in Figs. 23—25. It should be noted 
that the flow patterns upstream of the airfoil in Figs. 23 and 25 are 
significantly different with the flow for ot < 0 (Fig. 25) showing 
considerably more turning in the region upstream of the leading edge region. 
Further discussion of this case is presented in Ref. 36. 


NACA 4412 Airfoil 


The next calculation to be discussed is flow about an NACA 4412 
airfoil at high incidence. The calculation models the experiment of Coles 
and Wadcock (Ref. 65) which considered flow about an NACA 4412 airfoil at 
Reynolds number based upon chord of 1.5 x 10 and Mach number of 
approximately .07. The airfoil was placed in a wind tunnel and oriented at 
geometric incidence of 13.7 degrees to the free stream. Data taken included 
static pressure measurements on the airfoil surface as well as velocity 
profile measurements at selected boundary layer and wake locations. Details 
of the experiment can be found in the cited AIAA Journal article (Ref. 65). 
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The NACA 4412 calculation differed from those previously presented in 
two ways. First, the calculation was run with slip tunnel wall boundary 
conditions as discussed in the boundary condition section of this report and 
the calculation was run in a 141 x 39 mesh as opposed to the 81 x 39 mesh 
previously used. Finally, the calculation was run using an artificial 
dissipation parameter, a, of 0.05. As shall be shown, the increased 
resolution and decreased dissipation parameter had a major beneficial effect 
on the quantitative comparison between calculation and measurement. 

The NACA 4412 calculation was carried out with a grid containing 141 
pseudo-azimuthal (streamwise) grid points and 39 pseudo-radial grid points 
with high streamwise resolution being obtained in the vicinity of the airfoil 
leading edge and high normal resolution being obtained in the airfoil 
boundary layer. The streamwise grid spacing in the leading edge region was 
approximately .002 chords and the radial grid spacing at the airfoil surface 
was approximately 10 -5 chords. 

Calculations were made for flow at four incidences; 7.7 degrees, 

10.8 degrees, 12.3 degrees and 13.7 degrees. The calculation at 7 . 7 ° was 
initiated from a uniform flow with a no— slip condition gradually applied by 
decreasing the velocity over several grid points in the vicinity of the 
airfoil surface. The calculation was made using the time-conditioning 
methods suggested in Refs. 66 and 67. It should be noted that the zero lift 
incidence for the NACA 4412 airfoil Is —4° (Ref. 68) and, therefore, the 
cases considered represent, airfoil flow fields at large incidence. Predicted 
surface pressure distributions for each case are shown in Fig. 26. Although 
to our knowledge no data are available for comparison at the lower incidence 
angles, the comparison between data and prediction is excellent for the 13.7 
degree case. 

A comparison also can be made of the predicted and measured separation 
points. At the two lower incidences only a small amount of separation was 
predicted as the boundary layer separation occurred at x/c = .96. However, 
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at 12.3° incidence the separation point moved upstream to x/c = .82 and at 
13.7° it progressed to x/c = .72. This final value is in good agreement with 
the Coles-Wadcock data which indicated separation at x/c = .70. Although 
detailed pressure distribution comparisons cannot be made at the lower 
incidences, lift coefficient comparisons can be made to assess the prediction 
procedure. A comparison between predicted and measured lift coefficient is 
shown in Fig. 27 where the data are taken from Abbott and von Doenhoff 
(Ref. 68) for Re c = 3.0 x 10 6 . Although the predictions and data apply to 
somewhat different Reynolds numbers, the results should be insensitive to 
Reynolds number at lower incidences and only become somewhat sensitive as 
stall is approached. As can be seen in Fig. 27, the comparison is very good. 

The final comparisons concern velocity profiles. The Coles-Wadcock data 
give velocity profiles for both the velocity component parallel to the line 
connecting the airfoil leading and trailing edge locations and the component 
normal to this line. The streamwise velocity profiles are compared in 
Fig. 28, and the normal velocity profiles are compared in Fig. 29. 

Considering first the streamwise velocity profiles, it should be noted that 
the velocity is normalized by U re f, the velocity at a specified tunnel 
location and, therefore, U/U re f does not necessarily approach unity at the 
edge of the profile. In all cases, however, the predicted and measured edge 
velocities were in good agreement. Of the profiles compared in Fig. 28, two 
are taken in the aft region suction surface boundary layer and two in the 
wake. Both boundary layer profiles (x/c = .642 and x/c = .908) are in a 
strong adverse pressure gradient region. The predicted velocity profiles are 
somewhat thicker than the measured profiles and appear somewhat more advanced 
toward separation. However, the comparison is reasonably good in a 
qualitative sense. Similarly the wake profiles are affected by the somewhat 
too large prediction of the suction surface boundary layer thickness, 
however, the wake development also is reasonably good. Obviously, further 
investigations concerning both grid resolution (particularly in the 
streamwise direction) and turbulence and transition modeling may be required 
if the comparison is to improve. In particular based upon some more recent 
studies for cascade flows accurate prediction of the transition point are 
required to obtain accurate predictions of the suction surface boundary layer 
profile. Comparisons between the predicted and measured normal velocity 
profiles are given in Fig. 29. Again, the predictions and measurements are 
in reasonable agreement. 


33 


A more comprehensive view of the flow field can be obtained via vector 
and contour plots. A velocity vector plot of the flow field is shown in 

30 and contours of U— velocity, W— velocity and pressure coefficients are 
presented in Figs. 31—33. The velocity vector plot clearly shows the flow 
turning as it passes about the airfoil and the boundary layer build-up on the 
suction surface. The development of the large W-velocity component in the 
leading edge region is shown in Fig. 32 and the pressure field is shown in 
Fig. 33. The leading edge region is shown in detail in Figs. 34-37 which 
considers the flow in the initial 5% chord region. 


NACA 4412 - Grid Resolution and Artificial Dissipation Study 

As part of the NACA 4412 study, an effort was undertaken to assess the 
effects of the artificial dissipation parameter, a, and grid resolution. As 
previously discussed, the present prediction procedure utilizes second order 
central spatial differences and consequently requires addition of artificial 
dissipation terms to suppress spatial oscillations. The present approach 
adds a second-order dissipative term in the following manner. During the 
calculation, the flow field is examined at each grid point and to maintain 
stability terms of the form 


^x and JLjfc 

P n 9 x 2 p n 3 y 2 

are added to the equation. The variable <f> is taken as p in the continuity 
equation, u in the x-momentum equation and v in the y-momentum equation and n 
is 1 for the continuity equation and zero for the momentum equations. The 
factor d x is taken as the maximum of opuAx-p and zero and the factor d y 
is taken as the maximum of opvAy— p and zero where o is the artificial 
dissipation parameter, p is density, u and v are velocity components, A x 
and A y are grid spacing and P is viscosity. These points are discussed in 
detail in the artificial dissipation subsection of the present report. 
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In addition to the calculations for a = .05, presented in the previous 
subsection, NACA 4412 airfoil calculations were also made for a = . 10 and 

0.5. The calculation for a = 0.10 was initiated from the o = .05 calculation 

by abruptly changing the value of the dissipation parameter as the solution 
was restarted. The prediction for a = .05 and .10 were nearly identical as 
shown in Fig. 38. The insensitivity of predicted flow to changes in a 
between 0.10 and 0.05 has since been observed in other calculations performed 
at SRA (e.g.; Refs. 62 and 69). After completing the a = . 10 calculation, 
the procedure was restarted with a = .50. This calculation proceeded to 
diverge from the o = .05, .10 results and the o = .50 calculation was 

discontinued after fifty time steps although it had not yet converged and was 

continuing to move farther from the previous solutions. Since these results 
clearly demonstrate the effort of increasing the artificial dissipation 
factor to 0.5, no attempt was made to obtain the converged solution. 

The results of Fig. 38 clearly indicate that for the case considered 
(i) once o = . 10 a further decrease to o = .05 does not affect the solution, 
and (ii) a value of a = .50 leads to an inaccurate solution and calculations 
run with this value of artificial dissipation must be regarded as suspect in 
terms of quantitative results. 

The second item considered under the present study concerns grid 
resolution. The NACA 4412 results presented in the previous subsection were 
obtained with a highly stretched grid having 141 pseudo-azimuthal grid points 
and 39 pseudo— radial grid points. The grid was stretched so as to obtain 
maximum radial resolution near the airfoil surface (the grid spacing here was 
10 -5 chords), and maximum streamwise resolution in the leading edge region 
(the grid spacing here being .002 chords). For the grid resolution study, 
the calculation at 7.7 degrees incidence was repeated for two new grid 
sizes. In each case, the pseudo— radial grid was kept the same (39 points) 
since it is felt that this grid is required due to the presence of the 
Heimenz layer in the leading edge region and the laminar sublayer on the mid- 
chord and trailing edge regions. However, the streamwise grid requirements 
are somewhat less definite and consequently the 7.7° incidence calculation 
was repeated for 81 and 51 streamwise grid points. In each case, the 
relative stretch was kept the same as used for the 141 grid point cases. The 
calculations for the three cases showed no significant differences. The 
surface pressures were nearly identical with some slight smearing of the 


35 



suction peak for the 81 and 51 point calculations# However, the predicted 
lift coefficients only varied by + 3 per cent and the skin friction 
coefficients were also in very close agreement. Therefore, for the 7.7° 
incidence case, a 51 x 39 grid appeared to be adequate. 

The 7.7° incidence case is not a severe test since the flow is 
reasonably well behaved in this case. Therefore, a second study was made for 
a case at 12.3° incidence. In this case, a converged surface pressure 
distribution was obtained with the 141 x 39 grid and the results are shown in 
Fig. 39. However, when the calculation was repeated for a 51 x 39 grid, the 
surface pressure distribution would not converge and a nonconverged result Is 
shown in Fig. 39. Obviously, in this more demanding case at higher 
incidence, 51 pseudo-azimuthal grid points are inadequate. 
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High Incidence Dynamic Stall Calculation - NACA 0012 Airfoil 


The final case considered is that of an NACA 0012 airfoil in dynamic 
stall. The calculation was run against the data of St. Hilaire, Carta, Fink 
and Jepson (Refs. 70 and 71). The case chosen was case 51.005 for a 
NACA 0012 airfoil oscillating in pitch. The free stream Reynolds number was 
2.08 x 10^, the free stream Mach number was 0.30, the reduced frequency was 
0.125, the mean incidence was 12°, and the oscillation amplitude was 8°. A 
converged calculation was first obtained for an airfoil at 4° incidence; this 
result is compared with the data of Gregory and O'Reilly (Ref. 72) in 
Fig. 40. As can be seen, the comparison for this steady calculation is very 
good. A comparison between predicted and measured lift and moment 
coefficients is presented in Figs. 41 and 42. In regard to the lift 
coefficient agreement is reasonably good on the upstroke, however, the 
measured lift stall appears to occur before the calculated lift stall. The 
dotted portion of the line represents that portion of the calculation for 
which it was necessary to increase the artificial dissipation factor to 0.5 
as will be discussed subsequently. During stall, the measured lift was less 
than that calculated. However, the agreement appears to be qualitatively 
reasonable for this case. Considering next the moment coefficient 
comparison, the agreement between calculation and measurement again is 
qualitatively good. 

Although instructive, the comparison of measured and calculated 
coefficients only gives a comparison of integrated quantities. 

Relatively small differences in pressure distributions can lead to 
significant differences in lift coefficient and large differences in moment 
coefficients. Obviously, more insight can be gained via a comparison of the 
pressure distributions. Such distributions are presented in Figs. 43-50. 
Three comparisons during the upstroke are shown in Figs. 43 and 44. As can 
be seen, the agreement is very good. The data were reconstructed from the 
Fourier coefficients given by St. Hilaire and Carta (Ref. 72). The third 
measured data point on the pressure surface (x/c « .066) gave very erratic 
results and was not plotted for most of the comparisons. The excellent 
comparisons shown in Figs. 43 and 44 give evidence to the time-accurate 
calculation for the surface pressure. 
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Figure 45 presents a comparison at a = 17.7°, ct > 0. This is near the 
incidence where stall would first be inferred fom the lift and moment curves 
of Figs. 41 and 42. The figure shows some discrepancy between predicted and 
measured values as the data present some evidence of a vortex being shed on 
the suction surface leading edge. The discrepancy increases in Fig. 46 where 
the data clearly indicate stall. The plateau in the calculation on the 
suction surface, x/c » .15, seems to indicate a vortex being initiated. 
Furthermore, the calculated maximum suction peak at a = 19.5°, Fig. 46, is 
considerably less than that at a = 17.7°, Fig. 45. Based upon the plateau 
and the drop in suction peak, the calculated distribution at 19.5° appears to 
be beginning the stall process. This agrees with the normal force and moment 
coefficients of Figs. 41 and 42. The data at 19.5° is presented with the 
calculation at a = 19.9°, a > 0 in Fig. 47. Although these are at different 
values of a, they represent pressure distributions at approximately the same 
incremental time after stall is initiated; distributions are remarkably 
similar. Figures 45-47 indicate that although the calculation predicts stall 
to occur after the measured incidence, once stall occurs the calculated and 
measured pressure distributions become quite similar. The major discrepancy 
in the calculated and measured values may lie in the prediction of vortex 
initiation. This in turn is influenced by turbulence and transition 
modelling. 

Comparisons over the downstroke are given in Figs. 48-50. Obviously, 
the basic trends are in agreement as a strong qualitative comparison is shown 
between the calculation and the measured data. Overall, the detailed 
pressure distributions show good agreement with data and present a more 
favorable comparison than would appear from the integrated coefficients 
presented in Figs. 41 and 42. Although surface pressure comparisons 
obviously do not represent the entire story, these results indicate 
substantial agreement between calculation and measurement for this very 
difficult test case. 
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Velocity vector plots and vorticity contour plots are presented in 
Figs. 51-54. The results are qualitatively similar to those obtained for 
previous cases. During much of the upstroke, the velocity field remains well 
behaved and the vorticity is confined to the immediate region of the airfoil 
surface. However, near the flow field changes dramatically as a large 

clockwise vortex appears on the suction surface. This is evident in both the 
vector plot. Fig. 51 and the vorticity plot Fig. 54. Also, note the 
appearance of the trailing edge counterclockwise vorticity. This phenomenon 
has been noted and discussed in detail by Robinson and Luttges in a study of 
airfoils in pitch (Ref. 74). By a = 18.3°, & < 0 the leading edge vortex has 
clearly broken away and the trailing edge vortex of opposite sign has 
increased in size and is moving somewhat upstream. This is demonstrated 
more clearly in the vorticity contour plot. By a = 16.5°, a < 0 the vortices 
are tending to interact and begin to move downstream. This process also has 
been discussed by Robinson and Luttges. During the initial portion of the 
downstroke the flow is a very complex one with interacting shed vortices of 
opposite sense. The grid used concentrates resolution near the airfoil; 
however, during this portion of the downstroke, the complex vortex 
interaction occurs in a region removed from the airfoil where grid resolution 
is relatively coarse. The calculation encountered stability problems which 
almost surely could have been solved by increased resolution in this region. 
However, generation of such a grid and the subsequent calculation was beyond 
the scope of the present effort and, therefore, the problem was addressed by 
increasing the value of a to 0.5 over a portion of the downstroke as shown 
in Figs. 41 and 42. Finally, by a = 9.5°, a < 0 the flow has fully 
recovered. 
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CONCLUSIONS 


A numerical procedure for calculating steady and unsteady airfoil flow 
fields has been developed, demonstrated and documented under the present 
contract effort. The procedure solves the full time-dependent, compressible 
Navier-Stokes equations via an alternating direction implicit (ADI) method. 
The calculations are performed on a highly stretched grid which places the 
first grid point off the airfoil in the viscous sublayer. No slip conditions 
applied at solid boundaries. The calculation has been run extensively 
with a mixing length turbulence model. Although not detailed in the present 
report, calculations have also been run with a two equation model for an 
NACA 0012 airfoil at 6° incidence. The procedure is capable of calculating 
steady solutions using a matrix preconditioning technique which allows rapid 
convergence over a Mach number range between virtually incompressible and 
transonic. Unsteady flows which require transient accuracy, can be made for 
M > 0.15. 

The procedure has been used for several calculations including steady 
flows about an NACA 0012 airfoil at zero and modest incidence and flow about 
an NACA 4412 airfoil at modest to high incidence. In all cases excellent 
agreement between calculated and measured pressure distributions was shown. 
High incidence unsteady calculations were made for an NACA 0012 airfoil at 
19° which showed good agreement with available data. Calculations of an 
airfoil in pitch below the stall angle also showed good agreement with data. 
Finally, a calculation was made for an airfoil in deep dynamic stall, 

4° _< cK 20°. This represents a very difficult case which contains dynamic 
effects, large scale separation and multiple interacting shed vortices. 
Considering the difficulty of the case, the agreement between calculated and 
measured pressure distribution was very good. 

In regard to boundary layer velocity profiles, the only comparison 
considered under this effort focused upon suction surface boundary layer on 
the aft section of a NACA 4412 airfoil. The calculated profiles agreed 
reasonably well with those measured, however, some discrepancies were 
apparent. Based upon other efforts at SRA the reason for the discrepancy 
appears to be due to the calculated transition location and prediction of 
transition location remains a subject of current investigation. 
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The calculation procedure gives very rapid convergence for cases in 
which a steady flow is sought. For the range of Mach numbers considered, 

.01 < Moo< .5 convergence is obtained within 70 time steps independent of 
Moo. It is anticipated that this convergence property would remain through 
the transonic region. The present code is a scalar code which requires 
15 secs/ time step for 5600 grid points on the CYBER 203. This translates 
into approximately 1000 secs/converged solution. Further code speed-up could 
reduce this run time by a factor of 10—20. For unsteady flow where 
transients must be resolved between 400 and 700 time steps are required per 
cycle. 

In summation, the code developed represents a powerful tool for 
calculating steady and unsteady airfoil flow fields. The calculations run to 
date indicate good agreement with pressure data even for the very demanding 
dynamic stall case. Agreement with velocity profile data although reasonable 
may require further efforts concentrating upon turbulence and transition 
modelling. 
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APPENDIX B - SOLUTION PROCEDURE [58] 

Background 

The solution procedure employs a consistently-split linearized block 
implicit (LBI) algorithm which has been discussed in detail in [29, 56]. 
There are two important elements of this method: 

(1) the use of a noniterative formal time linearization to 
produce a fully— coupled linear multidimensional scheme which is 
written in "block implicit" form; and 

(2) solution of this linearized coupled scheme using a consistent 
"splitting" (ADI scheme) patterned after the Douglas-Gunn [57] 
treatment of scalar ADI schemes. 

The method is thus referred to as a split linearized block implicit (LBI) 
scheme. The method has several attributes: 

(1) the noniterative linearization is efficient; 

(2) the fully-coupled linearized algorithm eliminates instabilities 
and/or extremely slow convergence rates often attributed to methods which 
employ ad hoc decoupling and linearization assumptions to identify 
nonlinear coefficients which are then treated by lag and update 
techniques; 

(3) the splitting or ADI technique produces an efficient algorithm 
which is stable for large time steps and also provides a means for 
convergence acceleration for further efficiency in computing steady 
solutions ; 

(4) intermediate steps of the splitting are consistent with 
the governing equations, and this means that the "physical" boundary 
conditions can be used for the intermediate solutions. Other splittings 
which are inconsistent can have several difficulties in satisfying 
physical boundary conditions [56]. 

(5) the convergence rate and overall efficiency of the algorithm are 
much less sensitive to mesh refinement and redistribution than algorithms 
based on explicit schemes or which employ ad hoc decoupling and 
linearization assumptions. This is important for accuracy and for 
computing turbulent flows with viscous sublayer resolution; and 
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(6) the method is general and is specifically designed for the 
complex systems of equations which govern multiscale viscous flow in 
complicated geometries. This same algorithm was later considered by Beam 
and Warming [75], but the ADI splitting was derived by approximate 
factorization instead of the Douglas— Gunn procedure. They refer to the 
algorithm as a "delta form" approximate factorization scheme. This scheme 
replaced an earlier non-delta form scheme [76], which has inconsistent 
intermediate steps. 


Spatial Differencing and Artificial Dissipation 


The spatial differencing procedures used are a straightforward 
adaption of those used in [29] and elsewhere. Three-point central 
difference formulas are used for spatial derivatives, including the 
first-derivative convection and pressure gradient terms. This has an 
advantage over one-sided formulas in flow calculations subject to 
"two point" boundary conditions (virtually all viscous or subsonic flows), 
in that all boundary conditions enter the algorithm implicitly. In 
practical flow calculations, artificial dissipation is usually needed and 
is added to control high-frequency numerical oscillations which otherwise 


occur with the central-difference formula. 

In the present investigation, artificial (anisotropic) dissipation 
terms of the form 2 

V 

( 1 ) 


E 

j 


d. 




h i 2 


3x-j 2 


are added to the right-hand side of each (k-th) component of the momentum 
equation, where for each coordinate direction x j , the artificial 
diffusivity dj is positive and is chosen as the larger of zero and the 
local quantity U e (a Re^ x -1)/Re. Here, the local cell Reynolds number 
Re^xj for the j-th direction is defined by 


Re Axj ‘ Re | PU j| V". (2) 

This treatment lowers the formal accuracy to 0 (Ax), but the functional 
form is such that accuracy in representing physical shear stresses in thin 
shear layers with small normal velocity is not seriously degraded. This 
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latter property follows from the anisotropic form of the dissipation and 
the combination of both small normal velocity and small grid spacing in 
thin shear layers. 

Split LBI Algorithm 
Linearization and Time Differencing 

The system of governing equations to be solved consists of three/four 
equations: continuity and two/ three components of momentum equation in 

three/four dependent variables: p, u, v, w. Using notation similar to 

that in [29], at a single grid point this system of equations can be 
written in the following form: 

3H( <f>)/3t = D($) + S( <J>) (3) 

where <(> is the column-vector of dependent variables, H and S are 
column— vector algebraic functions of <j>, and D is a column vector whose 
elements are the spatial differential operators which generate all spatial 
derivatives appearing in the governing equation associated with that 
element. 

The solution procedure is based on the following two-level implicit 
time-difference approximations of (3): 

(H n+1 - H n )/At = f5(D n+1 + S n+1 ) (1-3) (D n + S n ) (4) 

where, for example, H n+ 1 denotes H( <J» n+ 1 ) and At = t n+ l - t n . The 
parameter 3 (0.5 - 3 - 1) permits a variable time-centering of the scheme, 
with a truncation error of order [At^, (3 - 1/2) At]. 

A local time linearization (Taylor expansion about <|> n ) of requisite 
formal accuracy is introduced, and this serves to define a linear 
differential operator L (cf. [29]) such that 

D n+1 = D n + L n ($ n+1 - $ n ) + 0(At 2 ) (5) 

Similarly, 

H n+1 = H n + (3H/3<J>) n (<J> n+1 - ,j> n ) + 0 (At 2 ) (6) 

S n+1 = S n + ( 3 S/ 3 <(>) n (<|> n+1 - <j> n ) + 0 (At 2 ) (7) 

Eqs. (5-7) are inserted into Eq. (4) to obtain the following system which 

is linear in <J> n+ l 

(A - BAt L n ) ($ n+1 - <fr n ) = At (D n + S n ) (8) 
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and which is termed a linearized block implilcit (LBI) scheme. Here, A 
denotes a matrix defined by 

A = (3H/3<j>) n - 0At (3S/3<j>) n (9) 

Eq. (8) has 0 (At) accuracy unless H = <J>, in which case the accuracy is 
the same as Eq. (4). 

Special Treatment of Diffusive Terms 

The time differencing of diffusive terms is modified to accomodate 
cross - derivative terms and also turbulent viscosity and artificial 
dissipation coefficients which depend on the solution variables. Although 
formal linearization of the convection and pressure gradient terms and the 
resulting implicit coupling of variables is critical to the stability and 
rapid convergence of the algorithm, this does not appear to be important 
for the turbulent viscosity and artificial dissipation coefficients. 

Since the relationship between u e and dj and the mean flow variables 
is not conveniently linearized, these diffusive coefficients are evaluated 
explicitly at t n during each time step. Notationally , this is 
equivalent to neglecting terms proportional to 3 Pe/H or in 

L n , which are formally present in the Taylor expansion (5), but 
retaining all terms proportional to p e or dj in both L n and D n . 

It has been found through extensive experience that this has little 
if any effect on the performance of the algorithm. This treatment also 
has the added benefit that the turbulence model equations can be decoupled 
from the system of mean flow equations by an appropriate matrix 
partitioning [56] and solved separately in each step of the ADI solution 
procedure. This reduces the block size of the block tridiagonal systems 
which must be solved in each step and thus reduces the computational 
labor. 

In addition, the viscous terms in the present formulation include a 
number of spatial cross-derivative terms. Although it is possible to 
treat cross- derivative terms implicitly within the ADI treatment which 
follows, it is not at all convenient to do so; and consequently, all 
cross-derivative terms are evaluated explicitly at t n . For a scalar 
model equation representing combined convection and diffusion, it has been 
shown by Beam and Warming [77] that the explicit treatment of 
cross— derivative terms does not degrade the unconditional stability of the 

present algorithm. To preserve notational simplicity, it is understood 
that all cross-derivative terras appearing in L n are neglected but are 
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retained in D n . It is important to note that neglecting terms in L n 
has no effect on steady solutions of Eq. (8), since <j> n+ ^ - tp n = 0, 
and thus Eq. (8) reduces to the steady form of the equations: 

D n + S n = 0. Aside from stability considerations, the only effort 
of neglecting terms in L n is to introduce an 0 (At) truncation 
error. 

Consistent Splitting of the LBI Scheme 

To obtain an efficient algorithm, the linearized system (8) is split 
using ADI techniques. To obtain the split scheme, the multidimensional 
operator L is rewritten as the sum of three "one-dimensional" 
sub-operators (i = 1, 2, 3) each of which contains all terms having 

derivatives with respect to the i-th coordinate. The split form of 
Eq. (8) can be derived either as in [29, 56] by following the procedure 
described by Douglas and Gunn [57] in their generalization and unification 
of scalar ADI schemes, or using approximate factorization. For the 
present system of equations, the split algorithm is given by 

(A - PAtL”) (<t>* - <|> n ) = At (D n + S n ) (10a) 

n ** n * n 

(A - 6AtL 2 ) (<j> - <p ) = A ( - <j>) (10b) 

(A - BAtL^) (<}> n+1 - <J, n ) = A (<{>** - <j> n ) (10c) 

where <J> and <p are consistent intermediate solutions* If spatial 
derivatives appearing in Lj_ and D are replace by three-point difference 
formulas, as indicated previously, then each step in Eqs. (10a-c) can be 
solveby a block-tridiagonal elimination. 

Combining Eqs* (lOa-c) gives (11) 

(A - 3AtL“) A _1 (A - BAtL 2 ) A" 1 (A - BAtL^ ) (<|> n+1 - <}> n ) = At (D n + s") 

which approximates the unsplit scheme (8) to 0 (At 2 ). Since the 

intermediate steps are also consistent approximations for Eq. (8), 

physical boundary conditions can be used for <j>* and <{>** [29, 56]. 

Finally, since the are homogeneous operators, it follows from 

Eqs. (lOa-c) that steady solutions have the property that 

<j>n+l - <j>* = <{>** = <j)fi and satisfy n n 

D n + S n = 0 (12) 

The steady solution thus depends only on the spatial difference 

approximations used for (12), and does not depend on the solution 

algorithm itself. 
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Figure 1 - Sketch of constructive coordinate system. 
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Pressure coefficient 






T=3 . 08 

T=3 .78 

T=4 . 43 



Streamwise lqcation, x/c 


Figure 5 - Pressure distribution for 19° airfoil after cessation of 

airfoil motion (airfoil motion ceases at T=1.83). 



T=4 .43 

T=5 . 03 

T=5. 38 



Streamwise location, x/c 


Figure 6 - Pressure distribution for 19° airfoil after cessation of 

airfoil motion (airfoil motion ceases at T=1.83). 


O Data of Young, Meyers, and Hoad 

(NASA Technical Paper 1266 

AVRADCOM Technical Report 78-50) 

T=4 . 0 

T=4 .5 

T=5 . 0 

T=5. 6 

T-5.7 

T=5. 9 



Figure 8 - Experimentally measured velocity field, a = 19°. 

(Data of Young, Meyers and Hoad) 
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Figure 18 - Static pressure coefficient contours, a = 19° 


Me = ,5. 


Lift coefficient 


1.0 



Incidence, a 


Figure 19 - Lift vs. incidence curve for NACA 0012 airfoil in pitch, k - 0.25. 
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Surface pressure coefficient. 
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- Boundary layer and wake transverse velocity profiles for NACA 4412 airfoil 


Figure 29 








Figure 30 - Vector plot NACA 4412 airfoil, a = 13.7°. 
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Figure 36- W-Velocity contours, NACA 4412 airfoil, a = 13.7°. 
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Surface pressure coefficient. 




Pressure Coefficient 



Figure 39 - Grid resolution study. NACA 4412 airfoil 
at 12.3 degrees incidence. 
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Pressure coefficients 



Figure 40- Comparison of calculated and 


measured pressure distribution, NACA 0012, 0 - 4 ° 




Lift coefficient 



Figure 41 - Lift coefficient - NACA 0012 airfoil. 
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Moment coefficient 



0 4 8 12 16 20 7L‘ 


Incidence, a 


Figure 42 - Moment coefficient - NACA 0012 airfoil. 
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Pressure Coefficient, Cp Pressure Coefficient f Cp 


Calculated 




Figure 43 - Pressure coefficient comparison - NACA 0012 airfoil, 
a = 6°, 12°, d > 0 
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Pressure Coefficient, Cp 



Figure 45 - Pressure coefficient comparison - NACA 0012 airfoil, 
a = 17.7° , a > 0 
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Pressure Coefficient, Cp 



Figure 46 - Pressure coefficient comparison - NACA 0012 airfoil, 
a = 19.5°, a > 0 
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Pressur 





Pressure Coefficient, Cp 



Figure 48- Pressure coefficient comparison - NACA 0012 airfoil, 
a = 18.3°, a < 0 
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Figure 51 - Velocity vector plot - NACA 0012 airfoil. 
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a = 9-5° , a < 0 


Figure 52 - Velocity vector plot - NACA 0012 airfoil. 
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18.3°, a < 0 


city contours - NACA 0012 airfoil. 
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